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ABSTRACT 


iii 


To analyze the flow in a supercritical channel Junction, two dimensional, depth av- 
eraged unsteady flow equations are solved by using false transient approach. An 
algebraic coordinate transformation is used to obtain a rectangular computational 
domain to simplify the numerical treatment of boundaries. Results of numerical 
experimentation on 22.5° and 45° junctions are compared with those obtained by 
using approximate method of Hager as well as with the available experimental data 
. Comparison of the present numerical results for large confluence angles indicate 
satisfactory prediction of the maximum depth of flow when the flow in the main 
channel is dominant. However, for small confluence angles the results are not very 
satisfactory due to the inherent assumptions of the shallow water equations. The 
present numerical simulation is able to adequately simulate the water surface con- 


tours at the confluence qualitatively. 
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Chapter 1 


Introduction 

Open channel flows.whether in natural channels or in man-made situations, is pre- 
dominantly in subcritical mode. Supercritical flows are rather rare and are also 
more complicated to analyze than their subcritical counterparts. Not surprisingly, 
the knowledge base relating to these supercritical flow situations is relatively less. 

Channel junctions are necessary component of a channel network. Junction of 
flow of two or more supercritical flows are sometimes encountered in mountaneous 
stream networks and in road drainage works. The flow situations in the channel 
junctions may be any of the following kinds The flows in the lateral channel and 
main channel before the junction are supercritical but after the junction it becomes 
subcritical or flow in any of the channels is supercritical and after the junctions 
it becomes subcritical or it may be that it is supercritical in both the channels 
and r emains so after the junction. Here such type of junctions are being studied 
where the flows are supercritical in main channel, lateral channel and it remains 
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supercritical after the junction. Here onwards this type of junctions will be referred 
to as SC junctions for the rest of the report. Supercritical flow in a combining 
channel junction is accompanied by standing waves which persists for long distances 
downstream of the channel junction. The side walls for such a junction have to be 
considerably higher than normal because the maximum height of the free surface 
may be 10 to 20 times as large as the depth of incoming flow (Behlke and Pritchett 
1966, Greated 1968 ). A hydraulic jump may form for certain combinations of 
junction configuration and the inflow r conditions. This leads to the problem in closed 
channels such as tunnel spillways due to pulsation, air entrainment and transition 
to pressurized flow ( Hager 1989 ) 

Analysis of supercritical flow is more complicated than the subcritical flow be- 
cause of the presence of cross waves. Bowers (1950 ) was probably the first person 
to study SC junctions. Schnitter et al. ( 1955 ) conducted experimental studies in 
the channel junctions with small confluence angles and found that a significant cross 
wave pattern is developed when the Froude numbers in the two upstream branches 
are not equal. Behlke and Pritchett (1966) and Greated (1968) presented a ana- 
lytical approach for determining wave angles in simple junctions and identified wall 
pile up zones. Gerodetti (1978) conducted experiments in junctions of semi-circular 
channels and steep bottom slopes. Hager (1989) used both analytical and exper- 
imental approaches to analyze the main flow features in simple SC junctions. He 
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considered rectangular cross sections and junction angles of 22.5° and 45° . 

Characteristics of supercritical flow can also be obtained by solving the two 
dimensional steady flow equations numerically. In spite of the significant advances 
and the refinement of computations of open channel flows, supercritical flow still defy 

reliable analysis. The basic equations are hyperbolic partial differential equations 

J.,’" 

and^such method of charateristic may be used to solve these equations. Yet for the 
two dimensional cases the method is too complicated. Bagge and Herbich (1967), 
Herbich and Walsh (1972), Villagas (1976) and Dakshinamoorthy (1973, 1977) used 
this method for analyzing steady state supercritical flows in channel transitions. 
The analogy between gas dynamics and shallow water flow has prompted many 
researchers to borrow techniques developed in the former area and apply them to 
supercritical flows in open channels. There are two types of techniques to compute 
flow fields with shocks ; the shock fitting and shock capturing techniques. Shock 
fitting methods treat shocks as discontinuities and explicitly compute their motion 
and strength, while shock capturing techniques do not give any special treatments 
to the shocks and their motion and strength is predicted as the part of the solution. 

Jimenez (1987) and Jimenez and Chaudhary (1988) used shock capturing tech- 
nique to simulate supercritical flow in transitions. Shock capturing techniques are 
capable of solving steady supercritical flow equations and marching technique may 
be used to simulate the oblique jumps or shocks. It fails for the case when the flow 
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anywhere in the channel becomes sub critical as the governing equations become el- 
liptic. Using unsteady flow model to obtain steady state flow solution by solving 
false transient approach mixed subcritical and supercritical flow may r be simulated. 
Such techniques have been applied in the past to simulate supercritical flows in 
channel transitions, among others, by Ligget and Vasudev (1965) Pandoifi (1982). 
and Jimenez and Chaudhary (1988). Bhallamudi and Chaudhary (1992) used two 
dimensional unsteady flow equations and the false transient approach to solve the 
same problem. This model can also simulate mixed sub and supercritical flows 
satisfactorily. 

In the present study, a two dimensional numerical model based on the false tran- 
sient approach is presented for computation of supercritical flows in simple channel 
junctions. A simple algebraic technique is used to convert the physical domain of 
the channel junction into a convenient computational domain. This makes the ap- 
plication of shock capturing finite difference methods easier. The computed results 
are compared with the experimental and analytical results of Hager ( 1989 ) to 
demonstrate the validity of the model. Brief introduction of analytical approach is 


given in the next section . 



1.1 Analytical approach for maximum depth in 
the channel junction 

The following model is based essentially on the analytical approach of Hager (1989). 
Fig(l.l) shows an idealized channel junction which is considered to be divided into 

three zones, these are 

Lateral Channel Flow 



Fig(l.l): Supercritical flow in a channel junction 

1. Upstream zone ( suffix u ) 

2. Zone between angle 9 and ff ( suffix 1 ) 


3. Zone between wall and wave front with the angle {suffix 2 ) 
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Fig(l-l) also shows the streamline pattern of the main flow. The depth of flow 
increases across the jumps and reaches a maximum in the zone 2. This maximum 
flow depth and can be estimate, as - » 



Fig(1.2): Supercritical flow past a wall deflection 


follows. Consider a supercritical flow past a wall deflection as shown in the Fig(1.2), 
the following relationships are used to determine h 2 , F 2 , 8 2 if hi, Fi, u> are given 
(Subramanya.1995) 


y = 


h 

hi 



+ 8 Fi + sin 2 0 — 1 


h 2 _ tan/5 
hi tan(/5 - u) 

F?-±(y-i)(y + i) 2 



(1.1) 

( 1 . 2 ) 

(1.3) 


As the determination of h 2 , F 2 and 0 is implicit in the above system of equations, 
Hager and Bretz (1987) derived the following explicit approximate equation from 
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the above equations. 

h = V2F, sin J! - 1 

ti i 2 

8 = jj 4 j = — 

2%/2iq 

F 2 = cos 2 3 

2 v^sin/3(l - 5^) 

where 

/i = '/2F if 3 


(1-4) 

{1.5} 

( 1 . 6 ) 

(1-7) 


Applying above equations for region u and region 1 to get yi, iq. 3\ and then for 
region 1 to region 2 to get y 2 , Fo, 82, the depth ho can be computed. 

To calculate 6 Greated (1985) assumed uniform velocity and hydrostatic distri- 
bution and derived the following relationship. 


1 — y 2 = 2 y 2 F* sin 2 9 — 2 F? sin 2 (o! — 9) (1.8) 

Hager (1989) further simplified Eqn.(1.8) and expressed 9 as 



where y = ^ 

All the analysis by Hager is valid only if the confluence angle is < 45° 



Chapter 2 


Governing Equations 

Open channel flows are described by three conservation laws viz. laws of conservation 
of mass, momentum and energy. If there are no discontinuities like jump or bore 
the momentum coefficient ( 3 , and velocity coefficient a, axe same and hence the 
momentum equation and energy equation are identical (Cunje et al,1980). For a 
supercritical flow in a channel junction an oblique hydraulic jump takes place and 
hence the momentum equation is appropriate. Though shallow water equations are 
not exactly applicable to such a case yet more general equations describing flows 
in open channel like Bousinesq type equations (which include vertical acceleration 
effect ) are not simple even for one dimensional flows. Therefore, the shallow water 
equations are used in this study to analyze flows in SC junction. 

In this chapter the shallow water equations, alongwith their assumptions and 
limitations are presented. This is followed by the co-ordinate transformation proce- 
dure for transforming the physical plane into a rectangular plane. 
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2.1 Shallow Water Equations 

Although in reality the flow in a SC junction is three-dimensional, for the purpose 
of analysis a two-dimensional approximation is adopted. Towards this vertically 
averaged quantities are used which simplifies the analysis and yields the results to 
reasonable accuracy. 

The governing equations are derived by depth averaging the 3-dimensional equa- 
te mo item. 

tions^over the flow depth. The equations in cartesian system are written as 

Lt + F x + Gy + 5 = 0 (-•!) 


where 


U 


( h \ 

hu 

\ hv I 


F = 


ul i 

urh + \g\r 

uvh 


( 


G = 


\ 


\ 


/ 


vh 

urh 

u~h + kglr 
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/ 


\ 


gh(S ox Sj x j 

\ ~9^-{Soy ~ Sfy) J 


subscripts t, x and y in Eqn 2.1 represents the partial derivative with respect to 
time, longitudinal direction and lateral direction respectively. 

In the above equation, 
t = time 

u = velocity along longitudinal direction 
v = velocity along lateral direction 
h = water depth measured vertically 
g = acceleration due to gravity 
S ox = channel slope along longitudinal direction 
Soy = channel slope along lateral direction 
S fx = friction slope along longitudinal direction 

Sfy = friction slope along lateral direction 


The bottom slopes are estimated as: 


Sqx — sin oc x 


( 2 . 2 ) 
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Q x -,a y are the angles between the bottom of the channel and the longitudinal and 
lateral direction respectively as shown in Fig(2.1) and Fig<2.2). 

Friction slopes may be estimated using Manning's or Chezy’s equation. For 


example : 



where C = Chezy’s coefficient 


2.2 Assumptions: 


( 2 . 4 ) 

(2.5) 


Shallow water equations are based on the following assumptions : 

1. The pressure distribution in the vertical direction is assumed to be hydrostatic 
i.e. assuming the acceleration in the vertical direction negligible. 

2. Depth averaged values are used in two-dimensional modelling. 

3. Fluid is assumed to be incompressible and the density is same throughout. 

4. Velocity distribution is uniform along the depth. 

5. The channel bottom is rigid and the slope is negligible. 

6. Only shear stresses due to horizontal velocity components are significant, other 
viscous effects are neglected. 
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7. Frictional resistance is approximated by the Chezv's equation or Manning's 

equation. 

The chief assumptions viz. 1 and 2 are discussed in the following paragraphs: 

2.2.1 Hydrostatic Pressure distribution: 

The assumption of hydrostatic distribution is correct if the streamlines are parallel 
and straight. As long as high curvatures are not present the the pressure distribution 
is almost hydrostatic. Ligget (1975) had shown that the assumption is valid as long 
as a shallowness parameter h 0 /l 0 (h 0 is water depth and l 0 is characteristic length) 
is small. Later Jimenez (1987) concluded that the shallow water theory reasonably 
represents the supercritical flow if depth to width ratio is of the order of 0.1 and 
Froude number is not close to unity. The error is of the order of 20% if the depth 
to width ratio is increased to 0.2. This error is manifested in the wavelength of the 
resulting wave pattern. The hydrostatic assumption is not valid in the vicinity of 
the shock and some flow details are lost at the shock location. However the overall 
results are adequate for engineering purposes (Cunge 1975). 

2.2.2 Turbulent effects and Depth averaging: 

The depth integration of the three-dimensional equations of motion produce higher 
order terms known as effective stresses which are not included in the Eqn(2.l). Effec- 
tive stresses constitute laminar viscous stresses, turbulent stresses and stresses due 
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to depth averaging of advective terms. Laminar viscous stresses may be neglected 
but the other two are important. Consideration of turbulent stresses require the use 
of turbulence closure model which expresses stresses as a function e of the main flow 
variables. Rastogi and Rodi (1978), Puri and Kao (19841 and Leschziner and Rodi 
(1979) used K — e models for describing turbulence. Yreugdenhil and Wijbenga 
(1982) used constant eddy viscosity concept. However, the above closure models 
are tested only for subcritical flows and no information is available for supercriti- 
cal flows. Without considering turbulent terms the side wall boundary layer effects 
are omitted. Ippen and Harleman (1950) analyzed the case of oblique jump and 
concluded that reduction of velocity across the boundary layer reduces the vertical 
acceleration at the point of wall deflection. This makes the shock front steeper than 
in the case when the boundary layer in not included. 

Ippen and Harleman (1950) studied the effect of non-uniform velocity distribu- 
tion also. They concluded that the effect of velocity distribution may be considered 
negligible.The measurement of velocity across the depth on the front and back side 
of the jump indicated that the momentum correction factor ranged from 1.007 at 
F r = 6.3 to 1.015 at F r = 2.0. It should be noted that though the effect of magni- 
tude variation of velocity may not be significant yet the variation of the direction of 
resultant velocity across the depth is significant. Especially in the channel junction 
as found by Hager (1989) the velocity vectors at different depths for supercritical 
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flow are not same. Velocity is almost parallel to the side wall near the channel bot- 
tom and it is pointing towards the wall near the surface. Eqni’2.1) can not simulate 
this effect and the numerical results will be in error. 

2.3 Coordinate Transformation: 

In this study, a finite-difference method is used to solve Eqnf2.i l. Physical boundary 
as well as the finite difference grid in cartesian coordinate system for a channel 
junction are shown in Fig(2.3). 


Y 



Fig(2.3): Finte difference grid superposed on channel 

We notice that it becomes necessary to represent physical boundary in angular 
portion in a stair step fashion. It leaves sharp corners (marked 1 on left side and 2 on 
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right side ) which introduce numerical disturbances which spread quickly to entire 
mesh. This may result into numerical instability (Roache 1972 } particularly for 
supercritical flows which are hyperbolic in nature. An alternate method of arranging 
the grid could be as in Fig(2.4) 



Y 

X 


Fig (2. 4): A second arrangement of grid 


This again leaves wedges at junctions of two channels and also require different 
coordinate directions for/ straight and angular channels. These problems can be 
avoided by using a coordinate system such that the coordinate axes coincide with 
the boundaries (Anderson et. al. 1984 ). Use of a simple algebraic coordinate 
transformation suffices in the present case Fig(2.5a). The following transformation of 
independent variables convert the original physical plane Fig(2.5a) into a rectangular 


Fig(2.5a): Physical Plane 


Fig(2.5b): Computational Plane 
With this, coordinate axes coincide with physical boundaries (see Fig(2.5c)) 



l 


Fig(2.5c): Final Grid 
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The governing equations must be expressed in terms of the new independent 
variables for their application. This is done by applying the chain rule of partial 


differentiation as follows: 

d_ _ d_ d_ 

dx ~^ x d£ + Vx di] 
d 8 d 

dy~^ y d^ +Vy d?i 

Subscripts denote partial derivatives. Eqn(2.6) and Eqn(2.7) give 


( 2 . 8 ) 

(2.9) 




i; 


Vx = 0 

(2.10) 


= cot 

a; 


1 

Vy — * 

(2.11) 





sin a 


d 

d 

d 

d 

/ 1 \ d 

(2.12) 




= cot a— 

+ 1 7T- 

dx 


dy 

d£ 

Vsin aj or] 


Applying the aforesaid transformations to the Eqn(2.1) and subsequently rearrang- 
ing them yields the following equations. 


U t + F ( + Gr, + S = 0 (2.13) 

where 


U = 


hu 

J 


F = 


uh + hv cot « 
u 2 h + \gh 2 + huv cot a 
^ uvh + ( hv 2 + \gh 2 ) cot a j 
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G 


' vhj sin a 
uvh j sin a 
^ (v~h -f \gti~) / sin a J 


( 0 ^ 
gh{S ox Sf x ) 

y ghl.Soy Sfy) J 


Eqn(2.13) is solved in this study using a finite-difference method. The details of 
the numerical scheme are presented in the following chapter, 



Chapter 3 


Numerical scheme 

3.1 Introduction 

The shallow water equations are hyperbolic partial differential equations. Analytical 
solution for these equations is possible only for very simplified cases. Therefore one 
has to opt for the numerical techniques for solving these equations. Among numer- 
ical schemes either it may be finite-difference technique or finite-element technique. 
These days finite volume methods are also becoming popular to solve such type of 
problems yet these methods are used to solve the equations in the integral form . 
In this study a finite difference method is used to solve the Eqn(2.13) for a channel 
junction. It is required that we get the steady state conditions in the junction. For 
that we assume some initial conditions as constant depth and velocity throughout 
the main channel and constant velocity and constant depth in the lateral chan- 
nel. The boundary conditions are set equal to the steady state conditions. Then the 
two dimensional unsteady flow equations are solved for sufficient length of time until 



tNT9*L 

' rUR 




variation of flow conditions is negligible and the conditions have converged to stead 


state values for the given boundary conditions . The Method of characteristics with 
shock capturing procedure (for two dimensionals cases ! take high computational ef- 
fort. As such it is advantageous to use shock capturing techniques rather than the 
method of characteristics. In these methods the following three points should be 


noted.: 


1. The governing partial differential equations should be in conservation form 
because Eqn(2.1) expresses the law of conservation of mass and momentum, 
and these quantities are conserved through the jump. 


2. The numerical scheme should be conservative i.e it should maintain the dis- 
cretised version of conservation laws throughout the computational region. 


3. The numerical scheme should be dissipative in nature so that it could represent 
the energy loss associated with the hydraulic shocks. The dissipation should 
be in such a way that it is small in gradual variation, but comparitively large 
near the discontinuities (Abbot 1975). First order schemes are dissipative in 
nature but the higher order schemes which do not have this mechanism produce 
spurious oscillations near the discontinuities. These spurious oscillations are 
originated because the conservative variables being differentiated across the 
discontinuities (shocks) may have discontinuous derivatives. To reduce these 
oscillations some researchers generally go for artificial viscosity and align one 



of the coordinate axis in the direction of shock. 

Artificial viscosity which is added in order to reduce oscillations acts otherway and 
it adds to the dissipation. As the first order schemes are highly dissipative and the 
second order schemes are having dispersive error as well as these schemes produce 
spurious oscillations. Therefore need was felt to develop a scheme which is not highly 
dissipative as well as does not produce spurious oscillation. For that researchers 
thought of going for total variation diminishing (TYD) schemes. Although the 
approximate Riemann solver based on the flux difference splitting of Roe and a flux- 
vector splitting of Van-Leer approaches are very accurate yet the methods require 
significant computational and programming time because of the field by field 

1 



Fig ("3.1): Finite Difference Grid 

decomposition. In the present study non oscillatory numerical scheme based upon 
the Lax -Friedrich’s solver proposed by Marinko Nujic (1995) has been used. The 


scheme has an advantage of its simplicity and ease of implimentation. Further it the 
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scheme treats the dispersion equation or turbulence model equations in the same way 
because Lax - Friedrichs solver does not depend on the structure of the hyperbolic 
system like the approximate solver which use field by field decomposition. The finite 
difference grid used in the scheme is shown in Fig (3.1). Here Ac and A rj are the 
mesh steps in £ and n directions and i and j subscripts are used in the following 
section to indicate the position along £ and rj directions respectively. 



Fig(3.2): Finite Difference Grid with imaginary points 

The dependent variables are known at any time level either as initial conditions 
or through computations for the previous step. The values of the variables at time 
t + At are determined using the finite difference technique as detailed in the next 


section. 


3.2 Nujic’s Scheme 


The scheme presented by Xujic to solve the hyperbolic partial differential equation 
is second order accurate. Only the brief introduction of the scheme for this problem 
is being described here. For the equation (2.1 } the variable u may be approximated 
at t + At (i.e. n + 1 th time level) when everything is known at time t she. n' k time 
level) in the following predictor and corrector form. 


u * =u n _ f n \_^( g n g n A-S^At (3.1) 

ij Ax A i+ 2- J Jt ~ 3 ' Ay J 


At 

Ay 


= 0.5 


u7i + U* 


At 


(/* + i j - Hij) - fe+i 9 h-d s ^ At 


^ *- J Ax 


Ay 


(3-2) 


where / and g are fluxes in f and rj direction and f i+ ij and g i j+ 1 are defined as 


Ji + iJ = 0.5 [/5 + /y - 7 «j - «y)] 

= 05 [ ft .® + fty - 1 (“® - “0 1 


(3.3) 

(3.4) 


in which 


uf, = Uij + 0.5 5uij 


u«j = U {+ 1 j + 0.5dUi+i j 


(3-5) 

(3.6) 


where 

= minmod (u i+ ij — — wt-ij) 

<5u i+ ij = minmod (tq+ij - - “i+ij) 


(3.7) 

(3.8) 
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and 


ufj = Uij + O.odUij 
B ~ 

^ i.j 0 

where the quantities du id and dUi.j+i are 

5uij = minmod (u^i - u,^. u V-i - 

dtijj+i = minmod (ui,j+i - u { j, Ui, j+2 - Ujj+i) 

and the function minmod is defined as 


1 O. 


\a\ 

c ( 

!&! 

< ! 


3.10 


( 3 - 11 ) 

( 3 . 12 ) 


minmod{a. ft) = ^ b if jftj < jal and ab > 0 


0 i/ aft < 0 

v 

/jj , ftj , g^j and gfj can be calculated in the same way. 

The above procedure can be used for the interior nodes. Boundary- conditions 
have to be applied at the boundary nodes. 


3.3 Initial and Boundary conditions : 

In the false transient approach any initial conditions such as constant velocity and 
constant depth equal to the boundary value may be assumed for all nodes. The 
determination of u, v, h during the unsteady computations involve the application of 
boundary conditions. Hyperbolic equations are sensitive to the boundary conditions 
because errors introduced at the boundaries propagate and reflect throughout the 
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grid and in many cases the instability may result Anderson et. al. 1954 . There are 
two type of following boundaries in the present study . These are flow boundaries 
and solid side wall boundaries. 

3.3.1 Flow Boundaries 

Mathematical basis for the specification of the boundaries in the shallow water wave 
applications have been discussed by Stokes (1957) .Abbott ( 1979 ) and Cunje et. al. 
Since the flow is two dimensional supercritical flow three boundary conditions are 
to be specified at the inflow boundaries and none at the outflow boundary. 


Y 

A 



Fig(3.3): Illustration of Boundary Grid 

In the present study u, v, h are specified at the inflow boundaries of the main 
channel and lateral channel and the downstream nodes these quantities are obtained 

4 

by extrapolation from the valmyat the interior nodes already calculated . Since the 
flow is supercritical their extrapolation affects only a few points at the downstream 
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end and the disturbance can not travel upstream . 


3.3.2 Solid Side wall boundary 


Since all the stresses other than the bottom stresses are neglected the slip condition 
is the proper boundary condition for the solid boundaries . The basic requirement 
is no flow normal to the boundaries, which is expressed as 

tan 5 = — (3.13) 

u 

where S is the angle between the wall and the axis. For the side walls of the main 
channel reflection procedure has been adopted because v = 0 normal to the side walls 
as the walls are aligned along the x-axis. Because of the coordinate transformation 
the reflection procedure is slightly modified. The mathematical representations of 
reflection procedure in cartesian coordinate system is as below: 


dh 

dl/wall 

du 


= 0 


= 0 


d Uwall 
dv 2 v wa u 


(3-14) 

(3.15) 

(3.1G) 


dy wa ii Ay 

Substituting Eqn(2.6) to Eqn(2.12) in Eqn(3.14) to Eqn(3.16) results in the following 


dh dh _ 
dS, + dr] 

(3.17) 

du du _ 

OS, + di] 

(3.18) 


equations 
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Ov dv 2v 

C ° Sa ^ + ^ = A^ 


13 . 19 ) 




x - direction 


Fig(3.4) : Illustration of boundary Conditions 
at the wall in lateral channel. 

Finite difference approximation of equations 3.12 to 3.14 for the solid wall shown in 
Fij>(3.3) gives t he following explicit, equations for the determining h, u and v at any 

node, (i,l) 


A || 

Ili t l = /lj,2 + (/h+1,2 — hi. 2) COS Ot 

A?/ . v 

Uij = «i,2 + — — («*+l ,2 “ »i,2) ( OSfV 

A?? , v 

ViJ = —Vi. 2 + (^1+1,2 u i.2j COSO 


(3.20) 

(3.21) 

(3.22) 


As u, u ,/i and fluxes are required at wail which is at midway of the nodal point as 
shown in the Fig(3.3).so the depth and velocity is the average of the depth at the 
ficticious points and the nodes adjacent to walls and velocity v = 0. We get the 
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quantities at the wall as 


\ 71 

hvaii = hi£ 4- 0.5—: - ft- 2 ) cos a 

_ _ A i] , 

^ wall — Ut,2 T 0*0 “ {'Uj'4.12. — Uj.2) COSO 


( 3 . 23 ) 

( 3 . 24 ; 


For the side wall the condition c = 0 can not be applied . The velocity normal to the 
face should be zero , therefore on the side walls c: the lateral channel the boundary 
conditions suggested by Bhallamudi(1992) are considered as 


h 


■ ficticious 


= hi 


'■interior 


V fi 


ficticious 


= \ 


interior 


(3.25) 


(3.26) 


which means that the depth and magnitude of the resultant velocity at the ficticious 
points and the interior points are equal and the directions of the velocities are in 
such a way that it cancels the normal components . Finally we get the depth and 
velocities at the wall as: 

V = \f “interior + l&arwr ( 3 * 27 ) 

u wa ii = V cos (a — 6) cos a (3.28) 

Vwaii — V cos (a - 0) sin a (3.29) 


3.4 Stability Condition 

Courant-Friedrich-Levy condition has mathematical foundation for the linear equa- 
tion. Shallow water equations are nonlinear and hence this condition can be used 
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as a guideline for the nonlinear equations but can not give full justification for such 
equation. Besides that it is not fully justified for TYD schemes . The condition for 
two dimensional flows is given as: 


CN = 



At 


v / At 2 + A ?/ 2 
A£A?7 


< 1 


( 3 . 30 ) 


In terms of At 


At = 


CN 1 

v '+^\/iFT? 


( 3 . 31 ) 


where V is the velocity at the grid point and CN should be less than one. The 
scheme used in the present study is stable under some CFL conditions. In addition 
to CFL condition another factor 7 is also important for the stability of this scheme. 
For stability 7 should be such that 


7 > max|Aij| (3.32) 

f — 1) imax 
J ~ 1 j j max 

where A it j is the wave velocity at some node i, j 

It appears that in some cases the value of 7 may be relaxed and the behaviour 
of 7 has been studied by Marinko Nujic (not yet published). The value of 7 used in 
the present study is the value specified in Eqn(3.32). 


Chapter 4 


Varification of the Model 


4.1 Introduction 


The results of the mathematical model presented in the Chapter 2 and 3 are com- 
pared with the available laboratory test data ( Hager , 1989 ). The numerical results 
are also compared with the approximate analytical solutions 



Branch - 2 


Fig(4.1): Plan View of Experimental Arrangement 



32 


of Hager. The experiments were conducted by Hager in a rectangular . horizontal 
channel of 50 cm width. Two confluence elements were inserted in the channel 
as shown in the Fig(4.1), making the width of the three branches equal to 99 mm 
each. Two junction angles a = 22.5° and 45° were considered in this study. The flow- 
depths were controlled by two vertical, movable gates w-hich were placed immediately 
upstream of the junction . The length of the downstream branch was 60 cm and 
after this it again expanded to the full channel width of 50 cm. The main test section 
was situated between the gates and the end of the downstream branch. Two types 
of tests were conducted. In the first type, the flow- depths h u and hi, the velocities 
V u and Vi, the angles of the main wave front 9 and the maximum flow 7 depths h max 
were measured in the majority of the runs. The Froude numbers covered the range 
2.8 < F r < 16 for a = 22.5° and 3.3 < F r < 8.3 for a = 45°. In the second 
type of tests the entire three-dimensional velocity field and the flow- surfaces w r ere 
recorded for three selected runs. These selected runs were made for a = 22.5°. A 
point gauge (±0.1 mm ) and a propeller velocity meter (±3 cm/s ) were used to 
take measurements of depths and velocity respectively. Details of the experiments 


are available in Ref ( Hager, 1989 ). 


O O <j unpt>- 
O rb p p £2 zz 


Fig(4.2)b: Water Surface Contours for Run 1 (Observed) 



Fig(4.4)a: Water Surface Contours for Run 3 (Computed) 



Fig(4.4)b: Water Surface Contours 


for Run 3 (Observed) 


10*0L ^cnm 


36 


4.2 Study of Water Surface contours 


The entire water surface profiles were measured by Hager for the following sets of 

inflow conditions 


Run j 

Fu 

Fl | 

h u {cm) 

hi(cm) j 

1 

5.25 

6.2 

2.65 ! 

1.98 S 

2 

4.5 

4.5 

2.65 

2.65 | 

3 

6.2 

5.25 

2.65 

3.58 ! 


Table 4.1 : Details of Hager’s experiments of type 2 

The numerical model presented in the previous chapter was run for the above 
test conditions using a numerical grid of A£ = A 77 — 1.724 cm . This corresponds 
to 15 grids in the main channel. The computational step, At was controlled by the 
Courant number and grid size. The upper limit for the Courant Number to achieve 
n um erical stability required some amount of numerical experimentation and a value 
of 0.6 was adopted for these runs. In addition to the Courant number there is 
another parameter 7 in the Nujic’s scheme which controls the spurious oscillations 
and hence stability. The value of 7 for the above runs was taken to be equal to 
the maximum wave velocity . The roughness coefficient was assumed to be zero 
since no information on this was available in the literature. The computed water 
level contours for Runs 1,2 and 3 are shown in Fig(4.2a), Fig(4.3a) and Fig(4.4a) 
respectively. The measured water level contours are shown in Fig(4.2b), Fig(4.3b) 
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and Fig(4.4b) respectively. The comparison of computed and measured water level 
contours indicate that the nature of the water surface contour is simulated fairly 
well by the present numerical model. However the simulation of the water surface 
profile along the walls is not very good. This is reflected in the prediction of the 
maximum flow depth. This point is discussed in more detail in the next section. It 
has also been obsereved that the water surface contours are simulated much better 
if the depths in the two channels are of the same order. For example in Run 2 the 
depths in the main channel and the lateral channel are the same (2.65 cm) and the 
depth contours for this case are the best out of the three cases studied. It is thus 
concluded that the water surface contours are simulates qualitatively well by the 
numerical procedure of this study. 

4.3 Prediction of Maximum Flow Depth 

Tables 4.2 and 4.3 present the maximum flow depth predicted by the Hager's ap- 
proximate method and the present numerical model, and also as measured in the 
experiments for different inflow conditions for oc = 22.5° and o; = 45 respecti\elv 
A comparison of numerical and analytical results for a = 22.5° is shown in Fig(4.5). 
It can be observed that agreement is good for some runs while it is not very sat- 
isfactory for some others. Both the numerical and analytical solutions are based 
upon shallow water flow assumptions. However, the approximate method assumes 
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a straight line shock front, line AC (in Fig(l.l)}. However, the results from the 
numerical model, indicate that the main wavefront at the angle 9 does not remain 
straight but becomes curved. This is also pointed out by Hager. From Tables 4.2 
and 4.3 it can be noted that the results are better when F u h u jFihi is of the order 
of unity or the flows in the main channel is stronger than the flow in the lateral 
channel. Also the results are better for lower Froude number (F r < 7 ) especially 
when the Froude number in the lateral channel is not very high. 














Table 4.3: Table showing maximum depths computed numerically (h n ). 
analytically (, h a ) and determined experimentally (h e ) 
(confluence angle = 45°) 
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h u (cm.) 

hi (cm.) 

F u | 

Fi 

h n (cm.) 

h a (cm.) \ 

h e (cm.) ; 

1.240 

1.280 

4.640 [ 

4.370 j 

8.248 

7.33 

9.20 | 

1.250 

1.290 

5.800 

5.450 

10.390 : 

9.96 

12.30 

1.280 

1.320 

6.720 

6.050 

11.932 

12.20 

15.60 

2.480 

2.500 

4.080 

3.920 

14.279 

12.30 

15.70 

2.460 

2.520 

4.480 

4.320 

15.848 

13.90 ‘ 

17.40 

1.830 

2.460 

4.960 

4.170 

14.210 

12.40 

15.10 

1.790 

2.530 

5.540 

4.580 

16.335 

14.20 

17.20 

1.220 

2.460 

6.300 

4.320 

12.405 

12.30 

14.70 

1.230 

2.470 

7.000 

4.820 

16.252 

14.20 

18.40 

2.430 

1.850 

3.790 

4.340 

11.779 

10.50 

12.90 

2.460 

1.850 

4.280 

4.810 

12.981 

12.50 

15.60 

2.470 

1.840 

4.710 

5.250 

14.122 

14.10 

18.40 

2.510 

1.300 

3.790 

| 5.100 

9.534 

10.20 

12.90 

2.480 

1.280 

4.420 

5.930 

10.931 

12.20 

15.00 

2.470 

1.300 

4.940 

| 

6.550 

12.288 

14.20 

| 18.80 

2.450 

.680 

3.300 

6.580 

6. ISO 

7.70 

8.60 

2.480 

.700 

3.890 

6.980 

6.639 

9.20 

10.00 

2.480 

.690 

4.560 

8.300 

7.546 

11.11 

13.20 

.700 

.750 

5.530 

5.160 

5.684 

5.30 

7.70 

.710 

.760 

8.110 

7.400 

8.493 

8.89 

13.00 

.460 

.500 

5.410 

4.740 

3.534 

3.32 

4.50 

.470 

.510 

6.700 

6.120 

4.647 

4.59 

6.20 

.490 

.530 

7.070 

6.530 

5.147 

5.16 

7.60 

.510 

.560 

7.150 

6.530 

5.449 

5.46 

8.50 
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a-'" 3 ! A*v t l 0 ^ 1> 


Fig(4.5) Comparison of Numerical and 

o 

Analytical Results for angle=22.5 
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Fig (4.7) Comparison of Numerical and 

O 

Experimental Results for angle=22.5 
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Numerical 


Fig(4.8) Comparison of Numerical and 
Experimental Results for angle=45 
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Fij-(1.7) con i part's the* numerical results with the* experimental results for a - 
22.5 . The matching of the numerical and experimental results is satisfactory only 
for the small values of maximum flow depth. But in the majority of cases the pre- 
dicted maximum flow depth is underestimated by the numerical model . The error 
in the numerical prediction increases with the increase in the maximum flow depth. 
This indicates, as expected (Jimenez, 1987), that the shallow water assumptions 
breaks down when the flow depth in the channel is large compared to the channel 


width. The assumption of uniform velocity across the flow depth at a point may 
be another reason for the unsatisfactory performance of the numerical model. For 
example the three dimensional nature of the velocity field for Run 2 is shown in 
■ Fig(4.9>. As can be seen, both the magnitude as well as the direction of velocity 
change across the depth. 

The velocity vectors are seen to be parallel to the side walls in the two inflowing 
sections except near the confluence corner. Further in the junction, the flow near 
the outside walls remains parallel to the boundary, both at elevation z — 1 cm and 
s = 2.5 cm above the bottom. The streamlines turn to follow the main wave front 
at an angle 6. This turning is less near the bottom and becomes more pronounced 
towards the free surface. At z = 5 cm , the flow along the wave front turns toward the 
outer wall and impinges on it, creating stagnation conditions. A two dimensional 
model , whether numerical or analytical, can not simulate this effect completely. 
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This error due to the assumption of uniform velocity cross a vertical could also be 
significant in addition to the error due to the assumption of hydrostatic pressure 
distribution, (Sec. 2.2.1) 

The comparison of the numerical and analytical results for q = 45° is shown 
in Fig(4.6). It is clear that results of numerical study are better than the and the 
approximate method of Hager when the flow in the main channel is dominant (i.e. 
Fu.hu/ Fihi > 1). This is especially so when the Froude number in the main channel 
is higher. While the numerical study predicts better than the Hager's method, 
still the maximum depth is underestimated in relation to the experimental results 
(Fig(4.7)). The Prediction of the maximum depth by the numerical method is much 
better when a = 45° than for the case of 22.5°. Thus it is concluded that the 
present numerical simulation gives better results than Hagers approximate method 
for predicting the maximum depth of flow in the junction with large confluence 
angle (a = 45°). However, for smaller confluence angle (a = 22.5°) the results of 
numerical simulation are of comparable accuracy with those of Hager s method. 

4.4 Numerical Sensitivity of the Model 

The present numerical model is sensitive to the parameters like Courant number 
(CN) and the coefficient 7 • For the purpose of prediction of water surface contours 
CN = 0.6 gave satisfactory convergence. While testing the effect of CN value on the 
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convergence value of CN < 0.7 gave satisfactory results. The factor 7 should be the 
maximum wave velocity for a time step. Hence the value of maximum wave velocity 
was used as the appropriate value for yin the present study. Though Nujici 1989 
) states that the above criteria for 7 can be relaxed in some cases^ In the absence 
of any criteria to decide the value, no relaxation was attempted. It is felt that the 
value of CN to be adopted in a study would depend upon the criteria adopted for 
the selection of 7. Numerical experimentations are needed in a specific problem to 
identify the best combination of CN and 7 



Chapter 5 


Summary and Conclusions 

A two dimensional mathematical model is presented for analyzing the flow in a 
supercritical channel junction. An algebraic coordinate transformations is used to 
obtain a boundary fitted orthogonal computational grid. This makes it easier to 
apply the side wall boundary conditions. The unsteady governing equations are 
used in a false transient approach to simulate the steady supercritical flow. 

The salient conclusions of present study are as follows. 

1. Regarding the maximum depth of the confluence the numerical prediction is 
nearer to the experimental results when 

(a) the flow in the main channel is dominant ( F u h u / Fihi > 1 and especially 
when F u > Fi) 

(b) the depths in the main channel and the branch channel are of the same 
order of magnitude 

(c) the depth to width ratio is small and 
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(d) the confluence angle is large 

2. Regarding the Water Surface Contours the numerical prediction of the water 
surface contours at the junction agree adequately with the experimental obser- 
vation qualitatively. However, exact quantitative predictions at the boundaries 
can not be achieved due to the inherent nature of modelling. 

3. The present numerical simulation gives better results than Hager's approxi- 
mate method for predicting the maximum depth of flow in the junction with 
large confluence angle (a ~ 45°). However, for smaller confluence angle 
(a « 22.5°) the results of numerical simulation are of comparable accuracy 


with those of Hager’s method. 



xvjljj. rjiX£jL\K^Fjd 


Abbott, M.B., (1975), Weak Solutions of the Equations of Open C'h 
”, in Unsteady Flow in Open Channels , Mahmood,K. and Yevjevich 
Resources Publications, Chapter 7, pp. 283-311. 


anne! Flow 
V. (eds... Water 


Abbott, M.B., (1979), Computational Hydraulics; Elements of the Theory of Free 
Surface Flow , Pitman Publishing Limited, London. 

Anderson, D.A., Tannehill, J.D. and Pletcher, R.H., (1984). Computational 
Fluid Mechanics and Heat transfer, McGraw-Hill, New York. 

Bhallamudi, S.M. and Cliaudhary, M.H., (1992), “Computation of flows in 

Open-Channel Transitions”, Journal of Hydraulic Research, IAHR Vol 30 No 1 pp 

77-93. 

Bhargava, R., (1991), “Computation of Supercritical Flow in Channel Junctions”, 
thesis submitted for the degree of Master of Technology, Indian Institute of Technolog}', 
Kanpur. 


Behlke, C.E. and Pritchett, H.D., (1966), “The design of Supercritical Channel 
Junctions”, Highway Res. Record, No. 123, Publication 1365, 17-35, Highway Reasearch 
Board, National Research Council, Washington, D.C. 

Bowers, C.E., (1950), “Studies of Open Channel Junctions”, Technical Paper No. 6, 

Series B, Part v, Hydraulic Model Studies for Whiting Field Naval Air Station, Universit 
of Minnesota, St. Anthony Falls Hydraulic Laboratory. 

Cui*je, J.A., (1975), “Rapidly Varying Flow in Power and Pumping Canals, in 
Unsteady Flow in Open Channels”, (Eds. Mahmood, K. and Yevjevich, V.), Water 
Resources Publications, Fort Collins. 

Chaudhary, M.H. (1993), Open-Channel Flow, Prentice hall, New Jersey. 

Dakshinamoorthy, S.,(1973), “Some Numerical Studies on Supercritical Flow Prob- 
lems” , thesis submitted for the degree of Doctor of Philosophy, Indian Institute of Tech- 
nology, Kanpur. 

Dakshinamoorthy, S.,(1977). “High Velocity Flow through Expansions " ,17"' Congress 
IAHR, Baden-Baden, Vol. 2, pp.373-381. 

Gerodetti, M. ,(1978), “Schussrinnen im Strassenbau” , Strasse und Verkehr. Zurich, 
64(7), pp. 271-275. (in German). 


Greated, C.A. ,(1968), Supercritical Flow Through Junctions” 
23(8), pp. 693-695. 


La Houille Blanche. 


Hager, W.H. ,(1989), “Supercritical Flow in Channel Junctions”. 
Eng., ASCE, Vol. 115, No. 5, pp. 595-616. 


Journal of Hydr. 


Herbich, J. B. and Walsh, P., (1972), “Supercritical Flow in Rectangular Expan- 
sions, Jour. Hydr. Div., Amer. Soc. Civ. Engrs., Vol. 98, No. 9, Sept., pp. 1691-1700. 

Ippen, A.T. and Dawson, J.H., (1951), “Design of Channel Contractions". Sym- 
posium on High-Velocity Flow in Open Channels, Trans. Amer. Soc Civ En^rs Vol 
116, pp. 326-346. ‘ ° ^ ' 

Ippen, A. T. and Harlman, D. R. F. , (1950), Studies on the Validity of the 
Hydraulic Analogy to Supersonic Flow Parts I and II, USAF Technical Report No. 5985 
May. 

Nujic, Marinko, (1995), “Efficient implementation of non oscillatory schemes for 
computation of free-surface flow” journal of Hydraulic Research, Vol. 33. No. 1 pp. 
101 - 111 . 

Jimenez, O.F., (1987), “Computation of Supercritical Flow in Open Channels”, 
thesis submitted for the degree of Master of Science, Washington state University. Pull- 
man. 


Jimenez, O.F. and Chaudhary, M.H., (1988), “Computation of Supercritical 
Free-Surface flows”, Jour, of Hydr. Eng., ASCE, Vol. 114, No. 4, pp. 377-395. 

Liggett, J.A. and Vasudev, S.U., (1965), “Slope and Friction Effects in Two- 
dimensional High Speed Flow”, ll i/l Int. Congress, IAHR, Leningrad, Vol. 1, Paper 
1.25. 

Pandolfi, M., (1972), “Numerical Experiments on Free Surface Water Motion with 
Bores” , 4th Int. Conference on Numerical Methods in Fluid Dynamics, Lecture Notes in 
Physic, Springer- Verlag, pp. 304-312. 

Puri, A. N. and Kuo, C. Y. (1985), “Numerical Modelling of Subcritical Open 
Channel Flow Using The k-e Turbulence Model and The Penalty Function Finite Element 
Technique,” Appl. Math. Modelling, vol. 9, No. 2, Apr., pp.82-88. 

Rastogi, A. K. and Rodi, W., (1978), “Predictions of Heat and Mass transfer in 
Open Channels ,“ Jour. Hydr. Engr., Amer. Soc. Civ. engrs., vol. 104, no. 3, March, 
pp. 397-419 

Roache, P.J., (1972), Computational Fluid Dynamics , Hermosa Publishers. 

Schnitter, G., Muller, R. Caprez, V. and Bisaz, E., (1955), “Modellversuche 



fuer Kraftworkbauten im Wallis”, ausgefuehrt an dor Hydraulischen Abtcihmg der Vor- 
suchsanstalt fuer Wassebau und Erdbau an der ETH, Wasser-und Energie-wirtschaft. 
47(5-7). (in German). 

Stokes, J.J..(1957), Water waves, Interscience Publishers, New York. 

Subramanya, K. (1995), Flow in Open Channels , Tata McGraw-Hill, New Delhi. 

Villagas, F., (1966),' “ Design of the Punchina Spillway”, Water Power and Dam 
Construction, Nov., pp. 32-34 

Vreugdenhil, C.B, and Wijbenga, J.H.A,(1982) "Computations of Flow Pat- 
terns in Rivers,” Jour. Hydr. Div., Amer. Soc. Civ.Engrs., Vo). 108. No.ll. pp. 
1296-1310. 



